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f— ^ ' Moving grids are of interest in the numerical solution of hydrodynamical problems and in nu- 

merical relativity. We show that conventional integration methods for the simple wave equation in 
^SJ ■ one and more than one dimension exhibit a number of instabilities on moving grids. We introduce 

two techniques, which we call causal reconnection and time- symmetric ADI, which together allow 
integration of the wave equation with absolute local stability in any number of dimensions on grids 
ry\ • that may move very much faster than the wave speed and that can even accelerate. These methods 

allow very long time-steps, are fully second-order accurate, and offer the computational efficiency of 
operator-splitting. 

We develop causal reconnection first in the one-dimensional case: we find that a conventional 
implicit integration scheme that is unconditionally stable as long as the speed of the grid is smaller 
than that of the waves nevertheless turns unstable whenever the grid speed increases beyond this 
value. We introduce a notion of local stability for difference equations with variable coefficients. We 
show that, by "reconnecting" the computational molecule at each time-step in such a way as to ensure 
that its members at different time-steps are within one another's causal domains, one eliminates the 
instability, even if the grid accelerates. This permits very long time-steps on rapidly moving grids. 
^-» , The method extends in a straightforward and efficient way to more than one dimension. 

However, in more than one dimension, it is very desirable to use operator-splitting techniques 
to reduce the computational demands of implicit methods, and we find that standard schemes 
for integrating the wave equation — Lees' First and Second Alternating Direction Implicit (ADI) 
methods — go unstable for quite small grid velocities. Lees' first method, which is only first-order 
accurate on a shifting grid, has mild but nevertheless significant instabilities. Lees' second method, 
which is second-order accurate, is very unstable. 

By adopting a systematic approach to the design of ADI schemes, we develop a new ADI method 
that cures the instability for all velocities in any direction up to the wave speed. This scheme is 
uniquely defined by a simple physical principle: the ADI difference equations should be invariant 
under time- inversion. (The wave equation itself and the full implicit difference equations satisfy this 
criterion, but neither of Lees' methods do.) This new time-symmetric ADI scheme is, as a bonus, 
second-order accurate. It is thus far more efficient than a full implicit scheme, just as stable, and 
just as accurate. 

By implementing causal reconnection of the computational molecules, we extend the time- 
symmetric ADI scheme to arrive at a scheme that is second order accurate, computationally ef- 
ficient and unconditionally locally stable for all grid speeds and long time-steps. We have tested the 
method by integrating the wave equation on a rotating grid, where it remains stable even when the 
grid speed at the edge is 15 times the wave speed. Because our methods are based on simple phys- 
ical principles, they should generalize in a straightforward way to many other hyperbolic systems. 
We discuss briefly their application to general relativity and their potential generalization to fluid 
dynamics. 
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I. INTRODUCTION. 

In the numerical study of wave phenomena it is often necessary to use a reference frame that is moving with respect 
to the medium in which the waves propagate. This could be the case, for example, when studying the waves generated 
by a moving source, where it may prove convenient to use a reference frame attached to this source. In some cases, 
one may even need to use a frame that moves faster than the waves themselves, as in the case of a supersonic flow. 
In general relativity, especially in black-hole problems, one may have to use a grid that shifts rapidly, even faster 
than light. All these problems arise in more than one spatial dimension, where computational efficiency may make 
stringent demands on the algorithm. It is a common experience to find that standard algorithms seem to go unstable 
in realistic problems. In this paper, by studying the simple wave equation, we show that the consistent application of 
two fundamental physical principles — causality and time-reversal-invariance — produces remarkably stable, efficient 



and accurate integration methods. These principles can easily be applied to more complex physical systems, where 
we would expect similar benefits. 

Our principal motivation for studying these techniques is the development of suitable algorithms for the numerical 
simulation of moving, interacting black holes. Relativists have long acknowledged the importance of using shifting 
grids in some problems, but to our knowledge there has been no systematic study of the effects of such shifts on the 
stability of numerical algorithms. In the next two paragraphs we develop this motivation. Readers not concerned 
with numerical relativity may skip these without loss of continuity. 

Let us consider the requirements that black-hole problems will make of our algorithms. Within the context of 
the 3 + 1 formalism of General Relativity ( jjj , H ) , it would seem to be desirable to develop methods on a quasi- 
rectangular 3-dimensional grid, so that no special coordinate features would prevent one from studying quite general 
problems. If we imagine a picture in which a black hole moves "through" such a grid, much the way a star would if 
it were interacting with another, then some requirements become clear: 

1. Grid points will move from outside to inside the horizon, but the grid as a whole should not be sucked in. 
This may require an inner boundary to the grid, say on a marginally trapped surface, and this boundary will 
have to move at faster than the speed of light. Grid points may cross this boundary and be forgotten, at least 
temporarily, but others will emerge on the other side of the boundary. 

2. Grid points that so emerge will then move from inside to outside the horizon as the hole passes over them; this 
will require grids that shift faster than light. This is inescapable unless one ties the grid to the hole as it moves. 

3. If two black holes begin in orbit around one another, then it may be desirable to adopt a grid that rotates with 
respect to infinity, in which the holes move relatively slowly at first. In such a grid one would expect that one 
could take long time-steps without losing accuracy, since not much happens initially. One therefore would like 
to be free of the Courant condition on time-steps, i.e. one wants to use implicit methods. 

4. Integrating the equations of general relativity on a grid with reasonable resolution will tax the capacity of the 
best available computers for some time to come. Full implicit schemes are very time-consuming in more than 
one spatial dimension, because they require the inversion of huge sparse matrices. Alternating-direction-implicit 
(ADI) schemes reduce this burden enormously by turning the integration into a succession of one-dimensional 
implicit integrations, so an ADI scheme that can cope with grid shifts is very desirable. 

In this paper we show that it should be possible to develop stable methods that satisfy the last three requirements 
above: ADI schemes that are absolutely stable and computationally efficient even on grids that shift at many times 
the speed of light. As a bonus, our ADI methods preserve the second-order accuracy of the full implicit equations. 
The first requirement, that of dealing with an inner boundary that moves faster than light, is closely related to these 
techniques and will be addressed elsewhere. 

Having these requirements in mind, we have studied the effects that the use of a moving reference frame has on 
the finite difference approximation to the simple wave equation, centering our attention particularly in the stability 
properties. The wave equation is the simplest system, so the instabilities we find in the standard ADI methods should 
certainly also be present when they are applied to more realistic physical systems. Of course, the wave equation is 
much simpler than other systems, so it is possible that methods that stabilize its integration will not extend to other 
systems. However, the principles that we find here are of such a fundamental physical nature that it seems certain 
that they should be applied wherever possible. Other kinds of instabilities may of course arise in complex systems, 
especially those directly due to nonlinearity, but we feel that moving-grid instabilities are likely to be cured by the 
methods we describe here. 

We shall conclude this introductory section by summarizing the approach and results of the following sections. 



In the second section we develop the mathematical framework of shifting grids. Then in Section III we study the 
one-dimensional wave equation. We find simple implicit finite difference schemes that are locally stable for any speed 
up to that of the waves, even when the grid is accelerating as well as moving. When formulated on a grid that is 
moving, and even accelerating, it is not immediately obvious how one defines stability: solutions of the differential 
equation do not have simple harmonic time-dependence in this frame. We find that a satisfactory criterion for local 
stability of these simple schemes is that no solutions of the difference equations should grow faster anywhere on the 
grid than local solutions of the differential equation. 

However, as soon as the reference frame moves faster than the wave speed, these schemes become highly unstable. 
We trace the origin of this instability to the fact that the computational molecules no longer represent in an adequate 
way the causal relationships between the grid points. We find that by modifying the molecules so that they link 



a given point on one time-slice with one on the next one that is within the first point's cone of characteristics (its 
forward "light cone"), one can restore stability. We discuss one algorithm for doing this in Appendix B. 

We call this causal reconnection. It is important to note that this has a minimal impact on the integration scheme: 
for implicit schemes, the matrix that must be solved for the solution at a given time-step is constructed only from the 
relations between grid-points at that time-step, while causal reconnection affects only the relations between points on 
different time-steps. Thus, it can be incorporated into the part of the algorithm that constructs the "inhomogeneous 
terms" that generate the right-hand-side of the implicit matrix equation. For the 1-dimensional wave equation, the 
extra work involved in seeking out causally related grid points can be significant, but it becomes a smaller proportion 
of the overhead in more than one dimension, and for complicated systems of equations, such as one has in general 
relativity, the overhead will be a negligible fraction of the total work per time-step. We have tested causal reconnection 
and found it to be stable even on grids moving at many times the speed of the waves. It is also insensitive to the 
acceleration of the grid. 



We then move our attention in Section IV to operator-splitting ADI methods []4|, which are computationally efficient 
ways of implementing implicit schemes in more than one dimension. We find it helpful to derive ADI schemes from a 
more systematic point of view than one usually finds in expositions of this technique. The goal is to add extra terms 
to a set of difference equations that (i) do not change its accuracy, but (ii) replace the large sparse non-tridiagonal 
matrix which has to be solved in implicit schemes with a matrix that is a simple product of tridiagonal matrices of 
the 1-D implicit form for each dimension, which are easy to solve. The extra terms are related to the "left-over" terms 
that appear as the difference between the true operator and its factored replacement acting on the data values on the 
final time-step. These final-time-step terms must be eliminated. They arc in effect replaced by similar terms from 
earlier time-steps, which replacement makes no difference when At — > 0, but which removes them from the matrix 
inversion and allows them to be included as part of the inhomogeneous terms in the matrix solution. Then the new 
equations will be a valid approximation to the differential equation but can be solved by a succession of (rapid) 1-D 
tridiagonal matrix solutions. 

When subjected to the same stability analysis as we devised for causal reconnection, the standard ADI methods 
show instabilities even when the reference frame moves very slowly. The instability is most marked in Lees' second 
method, in which the extra terms added in are of second order and therefore do not degrade the accuracy of the full 
implicit scheme. The instability is also present, albeit more weakly, in Lees' first scheme, which is only first-order 
accurate. 

We trace these instabilities to the fact that the extra terms added in either of the standard methods break the 
time-reversal invariance exhibited by the original differential equation and by the full implicit difference equations. 
Demanding that the extra terms be time- symmetric uniquely determines an ADI scheme that is essentially a hybrid of 
Lees' first and second methods. This time-symmetric ADI method turns out to be fully stable for all grid shifts up to 
the wave speed. Although not built in as a requirement, the new method also turns out to be second-order accurate. 

The method can then be extended to grid speeds larger than the wave speed by a direct generalization of the causal 
reconnection approach developed for the one-dimensional case. We demonstrate this by performing an integration on 
a rotating grid whose edge moves faster than the wave speed. 

In Appendix A, we derive the wave equation in the accelerating coordinate system using the efficient tensorial 
techniques of relativity. In Appendix B, we discuss one method of implementing causal reconnection. 

II. THE WAVE EQUATION ON A MOVING GRID. 

The wave equation is a good testing ground for any new algorithms for hyperbolic systems. The equations governing 
many wave systems can be reduced to the standard wave equation, and its cone of characteristics has the causal 
structure of space-time. We shall use it to test methods for integrating hyperbolic systems on moving grids. 

We consider the wave equation in an arbitrary number of spatial dimensions n , 
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written in a standard inertial coordinate system denoted by (t,^ 1 ). 

We are interested in finding a finite difference approximation to this equation using a grid of points that moves with 
an arbitrary non-uniform speed. Moreover, we will assume that the speed of each grid point can change with time. In 
order to represent this situation, we need to introduce a second coordinate system (t,x l ) that will be comoving with 
the grid. We introduce these coordinates in the continuous case by a transformation of the form 



x i = x % e)- (n.2) 

We have not changed the time coordinate, so we assume that the identification of surfaces of constant time does 
not change. This is thus not the usual Lorentz transformation of special relativity, so there is no reason for the form 
of the wave equation to remain invariant. This will have the implication that, in finite differences, the time interval 
between t = const slices will be constant, independent of position. For problems in general relativity, this is somewhat 
of a restriction, but we do not feel it is a serious one. If the causal relations are properly taken into account, then a 
spatial dependence in the lapse function ought not to change our physical conclusions. 

In Appendix A we show that the wave equation takes the following form in the new coordinates (t, x l ): 

dx 1 dx k c dx l dt dx 1 c 2 dt 2 

The following quantities derived from the coordinate transformation appear in the last equation: 

1 Br 1 

^:=--^-> (H.4) 

c at 

V- d ? d ? m *s 

^ := ^Sfei' (IL5) 

1=1 

^-E^ap ( IL6 ) 

g := det ( 9lJ ) , (II.7) 

I-:— -L{1(^) + ^[W(^-W]}. (n-8) 

Each of these quantities has a physical interpretation, which we now explain. Readers familiar with these ideas 
may skip to the next section. 

The shift vector fi % gives the relationship between the two coordinate systems on nearby surfaces of constant time.n 
Let the line of constant {£*} have coordinates {x\} at the lower hypersurface and {x l t+dt } at the upper hypersurface. 



From the definition of the shift vector in Equation II. 4, it is clear that 



x l (¥,t + dt) « x l (£ j ,t) - cf3 l dt. (II.9) 

As we illustrate in Figure U\ if one starts at any given point at time t, then by time t + dt the {x 1 } coordinates will 
have shifted by an amount equal to the shift vector times cdt relative to the {£*} (inertial) coordinates. The shift 
vector (3 l will in general be a function of both {x\} and t. 

We now introduce the spatial metric tensor gij, which we have defined in Equation [II. q . Its name comes from the 
fact that the distance dl between two points whose coordinates differ by d£ l in the original coordinates and by dx 1 
in the shifting coordinates is given by the Pythagorean Theorem: 

n n 

dl 2 = ^d£'d£' = J2 ffijdzW (11.10) 

i=i i,j=i 



That this gives Equation [1.5 for g^ is readily seen by substituting the following transformation from d£ l to dx 1 



into the first version of the Pythagorean Theorem: 



^ = E^- ( IL11 ) 

i—l 



This is just the standard definition of the shift vector in the 3+1 formalism of Numerical Relativity 
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FIG. 1. Shift vector f3\ 



The next tensor that appears in the general form of the wave equation is the inverse metric tensor given by 



Equation II. 6. This is the matrix inverse of the metric tensor, 



,9 l3 9jk 



= 51 



(11.12) 



as can easily be seen by substituting Equations II. 5 and II. 6 into the above. 

The final quantity we need is T l , a measure of the acceleration of the shifting coordinates with respect to the old 



ones, given by Equation 11.8. We will leave the full derivation of T* to Appendix A, but to illustrate our interpretation 
of it as an acceleration term, we shall explicitly evaluate it in the case where the new coordinates are obtained from 
the inertial ones by a simple shift independent of position. Then the shift vector (3 l is only a function of time, and 
the spatial metric gij is just the unit matrix: 



^ = 0, 

dx J 



9ij 



= 5 



13 



It is not difficult to see that in this case the coefficients T* reduce to: 

r = - - *t . 

c dt 



(11.13) 



(11.14) 



Since the shift vector gives the speed of the {x 1 } coordinates, the last expression implies that the T l coefficients are 
essentially the acceleration. 

Notice that if there is no acceleration, the only essential difference from the normal wave equation is the transport 
term f3 ■ V<j>, which arises as well in hydrodynamical problems. We will see that the local stability properties of the 
algorithms we study are determined mainly by (3 l , not T l , which is one reason we expect our analysis to have much 
wider applicability than just to problems involving the wave equation. 

Having derived the form of the wave equation in our new coordinates, we now establish a grid for formulating 
difference equations in these coordinates. By assumption, we take the time-interval At between successive surfaces of 
constant time to be uniform (independent of position) and constant (the same for any pair of surfaces). We take each 
grid point to have a fixed spatial coordinate position x l , and for convenience we take the spacing between grid points 
Ax 1 to be uniform in each coordinate direction. As seen in the inertial frame, the grid deforms itself as in Figure @. 
The corresponding picture in the ^-coordinate frame looks much more regular (Figure ra). 




FIG. 2. Grid in original coordinates, showing true distances. 
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FIG. 3. Grid in new coordinates. 



III. THE ONE DIMENSIONAL CASE. 



A. Finite difference approximation. 

The one-dimensional wave equation allows us to study shifting grids in a relatively simple fashion. The added 
complication of extra dimensions will be treated in the next section. 

In one spatial dimension, the metric, shift, and acceleration coefficients reduce to scalar functions: 



gn(x,t) = g(x,t), ' 
P l {x,t) = (3(x,t), 

r l (x,t) = r(x,t). 



(iili) 



Because the metric scales the squares of the coordinate distances (Equation II.10| ), it is convenient to define the 
linear scale function s(x,t) by 



s(x,t) := y/g{x, t) = 
so that the spatial proper distance is given by 

d£ = s(x, t) dx 



dx 



(111.2) 



(111.3) 
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FIG. 4. Computational molecule. 



Using this expression, Equation II. 2 becomes 
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For the finite difference approximation to this equation we employ the usual notation: 

4 := c(>(iAx,jAt). 
We define the first and second centered spatial differences as: 

5x4 ■= 4+1 ~<t>i-l, 

%4 :=4+i-^4 + 4-i 



(III.4) 



(III.5) 



(III.6) 



It is important to note that with the last definitions (S x ) ^ <5 2 . We can also define analogous differences for the 
time direction. 



We now write the finite difference approximation to the differential operators that appear in Equation [II. 4 using 
the computational molecule shown in Figure 0. We have: 
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Ett , 



where E u is the truncation error whose principal part is: 

,2 



Bk, 



(Aty 

12 



4 ±\3 



(dt4>) 



Similarly, for the mixed derivative in space and time we find: 

d x 4 +1 -^4~' 



(d x d t <t>)l = 



A Ax At 



with: 



E r 



C 



E x t , 



3 J.V 



(A X y (%d t <py.+ (AtY{d x dU) 



For the second derivative in the x direction we use an implicit approximation of the following form: 



ml = i 



%4 +1 



2 A3- 1 



K<t>. 



(Ax) 2 (Ax) 2 



+ (1-00 



(Ax) 2 



E x 



(III.7) 



(Illi 



(III.9) 



(111.10) 



(III.ll) 



where 9\ is an arbitrary parameter that gives the weight of the implicit terms. If 6\ = the approximation is 
explicit, while if #i = 1 all the weight is given to the initial and final time-steps of the molecule in the figure. Note 
that the last equation is symmetric in time. The error for this x derivative is: 



B _.-y*t {?M - ti ££ 



e> 2 Ati + 



Finally for the first derivative in x we take: 



(d x ^)l 



S X <P 



j+ l 



8 x <t>: 
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+ (1-02) 
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E x 



(111.12) 



(III. 13) 



where we have used again an implicit approximation with a different parameter 62 ■ The truncation error E x is: 



E x — — 



(As)' 
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(At) 2 



{d.dUt 



(III. 14) 



We can now write down a second order finite difference approximation to Equation III. 4 



i- -n <- 
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bj +1 - 24>l + <f>. 
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+ (1 - 02) 



5 X & 



where p is the 'Courant parameter' |p| given by: 
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(111.15) 



(111.16) 



The coefficients {s, (3, T} appearing in Equation [11.15 should be evaluated at the point (i,j) that corresponds to 
the center of the molecule. 

To arrive at the final form of the difference equation we multiplied it through by (At) . This means that the overall 
truncation error is now 



£111.15 = O 



(Ax) 2 (At) 2 



O 



(Atf 



(III. 17) 



Equation 111.15 is well studied in the particular case when (3 
that, 



= r = and s = 1 B. It is important to note 
the above finite difference approximation will be 
. Therefore the use of implicit 



because we use centered differences in the transport term 
implicit whenever the shift vector is different from zero, even when 9\ = 62 
approximations for the spatial derivatives does not add any extra numerical difficulty. 

We shall need to know how muc h num erical work is involved in using the implicit scheme. Suppose there are N 
spatial grid points. Then Equation [11.15 is to be solved for the N values {4>l , i = 1, . . . , N} at the final time-step. 
The equation for index i relates three such values, at points {i — 1, i, i + 1}. The system of equations therefore has 
the matrix form 



Q x 



d+i 



fiP, <t> 



p-^ 



(111.18) 



where Q x is a tridiagonal N x N matrix, and the inhomogeneous term / is constructed from field values at the first 
two time-steps. Solving a tridiagonal matrix involves O (N) operations. Since we also need O (N) operations for the 
solution of an explicit scheme, we see that the use of an implicit method in one dimension will increase the number of 
operations per time-step by at most a multiplier, independent of the number of grid points. Against this, the implicit 
scheme for certain choices of 9\ and 62 can, on a fixed grid, take much larger time-steps, limited only by accuracy 
considerations. In the next section, we shall show that this property of the implicit scheme in one dimension can, 
with suitable modifications, be extended to grids that shift essentially arbitrarily fast. 



B. Local stability: definition and analysis of the implicit scheme on a shifting grid. 



It is well known W] that the implicit approximation to the wave equation can be made unconditionally stable in 
the case when [3 = T = and s = 1 by using an implicit parameter 9\ > 1/2 . We are interested in studying 
under what conditions this property is preserved in the case of a shifting grid. The shifting grid introduces a major 
difficulty: the coefficients in the equation generally depend on both position and time. This complicates the definition 
of stability. 

This difficulty means that an analytic stability analysis must be local: we will actually only consider the stability of 



the difference equation obtained from Equation III. 15 by, at each point (x,t), taking the coefficients to be constant, 
with the values corresponding to that point. We feel that this is not a very restrictive assumption, since in practice 
instabilities usually appear as local phenomena , with the fastest growing modes having wavelengths comparable to 
the grid spacing. Moreover, if the coefficients in the difference equation are not practically constant over a few grid 
points, then we are probably not approximating the original differential equation adequately anyway. 

We will start then by considering the nature of the solutions of the differential equation in a very small region 
around the point (x, t) . As usual, we look for a solution of the form: 



<p(x,t) 



%ott ikx 



(111.19) 



Substituting this in Equation III. 4 gives the following "dispersion relation" for a: 
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^r+ikT 
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The general solution for a wavenumber k is 

<f>(x,t) - 
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Z_e l 



(111.20) 



(111.21) 



Clearly, if T ^ 0, then one of the independent solutions will grow with time, but the other one will decay because 
of what we shall call the analytic boundedness condition: 



\ e la + t e la - t \ 



1. 



(111.22) 



This does not mean that the system is physically unstable, but only that in an accelerating coordinate system (T =/= 0) 
the wave equation does not have purely sinusoidal solutions. One can understand this intuitively in the following way: 
Consider a sinusoidal solution in a static coordinate system. From the point of view of an accelerating observer, the 
frequency of this solution will be changing with time (he will be seeing more and more crests per unit time). This 
change in frequency will have important local effects. As a crest approaches our accelerated observer, he will see 
the wave function rising faster than an observer moving at the same speed, but not accelerating, would. Hence the 
appearance of locally growing modes in our analysis. Similarly, after a crest is reached, the accelerated observer will 
see the value of the wave function falling faster than a uniformly moving observer would. It is not difficult to see that 
the difference between the growth rate in the first case and the decay rate in the other, as seen by our two observers, 
will be the same. This is the origin of the analytic boundedness condition given above. 

The presence of the growing modes is crucial for our local stability analysis. Since the solutions of the differential 
equation can grow with time, we can not ask the solutions of the finite difference approximation not to do so. What 
we are entitled to ask is for the numerical solutions not to grow faster than the corresponding normal modes of the 
differential equation. Our stability criterion is, therefore: a difference equation is locally stable if every solution for a 
given wavenumber k is bounded in time by a solution of the differential equation for the same k. 

Bearing this in mind, we now proceed to an analogous analysis of the solutions of the finite difference scheme. We 
look for a local stability condition around the point (n, m) by making the substitution: 



WO 



m ikn Ax 



(111.23) 



Substituting the last expression in the finite difference approximation (Equation III. 15 ) we find a quadratic equation 
in ip of the form: 



A (iP) 2 + Bip + C = 0, 



(111.24) 



with coefficients given by: 
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ip sin (k Ax) 



2(l-0i)p 2 



cos (k Ax) — 1 



(3- -f(cAt)T 



- 1 



(111.25) 
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cos (k Ax) — 1 



C 



The two roots of this equation are: 



ip(l-6 2 )(cAt)T sin(fcAa;) 
1 



(111.26) 
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cos (k Ax) — 1 



ip sin (fc Ax) 
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(cAt)r 



■0=1 



-B ± (B 2 - 4AC) 
2A 
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and the general solution of the difference equation is 

ff=e^[2 + (^r+z_(^_n. 

It is not difficult to see that the coefficients A and C have the property: 

\A\ 2 = \C\ 2 - 29 2 p 2 /3 (cAt) T sin(fcAx) , 
which implies: 
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(111.27) 
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(111.30) 



(111.31) 



This contrasts with the differential case, Equation [III. 22] , where the product of the magni tudes of the two funda- 
mental solutions was 1. Since the ratio \C/A\ depends on the value of k in Equation [11.30 , there will always exist 
wavenumbers for which the product |^»+^»_| exceeds 1. This would seem to be undesirable from the point of view of 
stability, but we can eliminate it as a potential problem by setting from now on 



0. 



(111.32) 



This means that we will use an implicit approximation only for the second spatial derivatives {9\ ^ 0) , and not for 
the first spatial derivatives (02 = 0) . Since from now on we will have only one 9 parameter, we will change notation 
now and define 9 := 9\ . The solutions of the difference equation now satisfy: 



l^+V-l = 1. 
Next we introduce the amplification measure M : 



M := max 

k 



(W,h/M 2 ) , 



(111.33) 



(111.34) 



and analogously for the solutions of the differential equation. The amplification measure bounds the growth in the 
magnitude of any normal mode in one time-step. Our local stability condition is then equivalent to 



< M A na, 



(111.35) 



where M^ um and A^Anai are the amplification measures for the finite difference approximation and the differential 
equation respectively. 
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FIG. 5. Stability on a static grid. In the left-hand figure, we treat the explicit scheme, where we find, as expected, that 
instability sets in for Courant parameter p/s > 1. On the right, we see that a fully implicit scheme (8 = 0.5) is stable for all 
time-steps, again as expected. 



When all the parameters are free to take any value, Equation 111.35 is very complicated, and it is then difficult to 
find its consequences analytically. We shall therefore study this equation numerically, in order to find regions of the 
parameter space in which the finite difference scheme is stable. 

First let us consider the case of a static grid, = , T = 0. This case has, of course, been studied analytically [Q, 
and it is known that if 9 < 1/2 the generalized Courant stability condition is: 



V.s 



< 



1 



(1 - 29) ' 



(111.36) 



while if 6* > 1/2 the scheme is absolutely stable. In Figure |j| we show graphs of both the numerical amplification 
measure (solid line) and the one corresponding to the differential equation (dotted line). We have only plotted the 
functions for k Ax = 7r because this turns out to be the worst case. The first graph shows how for 9 — the scheme 
is stable for values of the Courant parameter p/s smaller than one. However, when this parameter takes values 
slightly larger than one, the numerical amplification measure begins to grow very fast. In the second case we see 
that for 9 — 1/2 the scheme is locally stable for all values of the Courant parameter, in agreement with the known 
stability condition given above. 

The next group of graphs (Figure g) shows the effect of a uniform shift. In both graphs we have assumed that there 
is no acceleration (T = 0) , and we have taken 9 = 1/2 in order to avoid any instability of the type seen in Figure ra. 
The first of these shows that the scheme remains locally stable for all values of the Courant paranreter, even when 
the grid shift speed s(3 is very close to 1. However, in the second graph we see that, as soon as s/3 becomes larger 
than 1, the scheme turns unstable for all values of p. In this last case there is no stable choice of time-step. This is 
a very important property: The finite difference scheme becomes unconditionally unstable whenever the shift is faster 
than the speed of the waves. 

Finally, in Figures ffl and H we consider the effects of an accelerating grid for the particular case when: 9 = 1/2 , 
s (3 — 1/2 , and s 2 T Ax = 1 . In Figure M we show the behavior of the amplification measure for the finite-difference 
equation and the differential equation for two different normal modes (two values of k). As we expect, the amplification 
measure corresponding to the differential equation, M^ na is no longer 1 . For the first graph we have k Ax = 1 and 
for the second k Ax — 2 . For the smaller wave number (larger wavelength) the amplification measures for the 
differential and finite-difference cases are relatively close to each other. As the wavenumber increases, the finite- 
difference amplification measure falls further below that of the differential equation, so that the finite-difference 
scheme remains stable (although less accurate). Figure || shows a surface plot of (l^Ana — Mn u111 ) in the region: 
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FIG. 6. Stability on a uniformly shifting grid. The figure on the left has a grid speed 0.9 times the wave speed. On the right 
the grid moves at 1.1 times the wave speed. In both cases we have set 9 = 1/2 and T = (no acceleration). 



p/s e (0,2) 



kAx € (0,7r), 



We clearly see how (MA na — M Num ) > in the whole region. Since k Ax = it corresponds to the smallest wavelength 
that can be represented on the grid (A = 2 Ax) , we find that the finite-difference scheme will be stable for all modes. 
We have searched through other values of T , and we have found that, although the details of the graphs change, 
the qualitative behavior is preserved. The acceleration parameter T thus seems to have no important effect on the 
stability of the scheme. 



In summary, our stability analysis shows that the finite difference scheme given by Equation III. 15 will be locally 
stable for all values of the Courant parameter p if the following conditions are satisfied: 



> 1/2, 



• 9 2 = , 

• \s/3\<l, 



irre- 



levant . 



(111.37) 



The limit on (3 is inconvenient in many problems, where it is desirable to have grids shifting faster than the wave 
speed. We turn now to a method for removing this restriction. 

C. Causal reconnection of the computational molecules. 

1. Causality problem. 

The causal structure of a grid shifting faster than the wave speed is particularly clearly illustrated in the original 
(£, t) coordinates. In Figure |9Jcz, we see how, for a very large shift, the individual grid points move faster than the 



12 



o 
< 



2 



k Ax = 1.00 
s/S = 0.50 



s 2 T Ax = 1.00 
i5 = 0.50 




l 
p/s 



o 
< 



2 



k Ax = 2.00 
s/S = 0.50 



s 2 T Ax = 1.00 
i5 = 0.50 



l 
p/s 



FIG. 7. Stability on an accelerating grid. For two different modes, the finite-difference amplification measure (solid curve) 
lies below that of the differential equation (dotted curve). This means the finite-difference scheme is stable, at least for these 
modes. 



i3 = 0.50 



s/S = 0.50 



s 2 T Ax = 1.00 



(M 



ANA 



M 



NUM 




-^ p/s 



p/s £ (0,2) 



k Ax e (0,77) 



FIG. 8. Stability on an accelerating grid. Here we show a surface plot of the difference between the analytical and numerical 
amplification measures. This difference is always positive, which means that the finite-difference scheme is stable in the whole 
region. 



13 



waves, that is, they move outside the light-cone.f] Since the differential equation propagates data along this cone, it 
seems plausible that the instability found in the previous section arises in the fact that the difference scheme attempts 
to determine the solution at points on the final time-step using data that are outside the past light-cone of these 
points. 

This suggests that we should not build the computational molecules from grid points with fixed index labels, 
but instead use those points that have the closest causal relationship (Figure ^6). We shall now proceed to show 
analytically how such a reconnection can stabilize the scheme. 

In order to build this causal molecule let us consider then an individual grid point at the last time level. We look 
for that grid point in the previous time level that is closest to it in the causal sense. Having found this point, we 
repeat the procedure to find the closest causally connected point in the first time level. In Appendix B we give a 
simple algorithm for finding these points in an integration of the wave equation. The algorithm adapts easily to 
other linear hyperbolic systems. We shall return in a later paper to its generalization to nonlinear equations, and in 
particular to the case of shocks in hydrodynamics. First we consider the general constraints on the time-step that 
causal reconnection imposes, and then address the issue of how much extra computational effort causal reconnection 
may involve. 

2. The causal reconnection condition. 

Since we permit the grid to move with an arbitrary non-uniform speed, there is no reason that these causally con- 
nected points should be in a straight line in either the original inertial reference frame (£, t) or in the moving reference 
frame (x, t). In the moving coordinate system (x, t) the relationship among these three points may generically look 
something like that shown in Figure |l0| . 

Accordingly, we introduce a new local coordinate system (a/, t') adapted to the three given points. In order to do 
this, it is convenient to introduce the interpolating second order polynomial that can be obtained from these three 
points: 

P(t) = A {t ~^ ] +B(t-t )+x , (111.38) 

where to is the time at the central point of the molecule, x to the position of that point, and: 



Xt +At - 2x to + X to -At \ R _ f xtp+At - X t -At 

(At) 2 )' '~\ 2 At 



A ;= ( ^t 0+ At-^t -r^ -At ] B . = | ^o+^-^o-^ ) . (ni.39) 



We define the new local coordinate system adapted to the causal molecule by: 

x 1 := (x-x to )~P(t) , ) 

\ (111.40) 

t' := t. J 

It can easily be seen that this new coordinate system (x', t') moves with respect to the old one (x, t) with a speed 
B at time t = to , and with a constant acceleration A . Since in general the value of the coefficients A and B will 
change from molecule to molecule, the above change of variables must be repeated for each molecule. We assume that 
this can be done in a smooth manner; this may not be possible in a nonlinear system, if the characteristics depend 
on the solution. 



In the primed coordinate system, the differential equation has exactly the same form as before (see Equation III. 4 ) , 
except for the following substitutions: 

(]^p + A{t - to) + B , r^r-4. (IH.41) 



2 From here on we will adopt the language of relativity and refer to the characteristic cone of the hyperbolic equation as the 
'light-cone'. 
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Light cone 




(a) Non-causal molecule. 



Light cone 




(b) Causal molecule. 



FIG. 9. Causal computational molecule. In both figures, the dashed lines represent the light-cone, and the thick solid lines 
show the computational molecule. Figure (a) shows the usual molecule that follows the motion of the grid points. Figure (6) 
shows the reconnected molecule, where we pay attention to the causal structure instead. 
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FIG. 10. Causally connected grid points. 



In the same way, the finite difference approximation will have the same form as before (Equation III.15| ), except for 
the substitutions: 



P 



P- 



B 



A 

72' 



(111.42) 



where the term with (t — to) has disappeared because in this case the coefficients should be evaluated at the center 
of the molecule where t = to ■ 



Since the original finite difference approximation could be made stable as long as Equation [11.37 was satisfied (and 
9 > 1/2), the analogous condition for a finite-difference scheme adapted to the new local coordinates takes the form: 



1 



P- 



B 



> 0. 



(111.43) 



We will say that the three given points form a proper causal molecule when the last condition is satisfied. In order 
to find when this happens, we will start by defining the effective numerical light-cone of the point Xt as the region 
between the lines: 



•'• i ( / ) := ■'•/„ r ( -/5 ± ij c (i- t ) + I Tc 2 {t - t ) 2 



This numerical light-cone will coincide with the exact light-cone when: 

dx dx dt 

We will also define the axis of the numerical light-cone as the line: 

1 r 

x a (f) := x to - (3c (t - t ) + - Tc (t- to)" 



(111.44) 



(111.45) 



(111.46) 



We will now show that if xt +At and Xt -At are inside the numerical light-cone of xt then the three points will 
form a proper causal molecule. From the definition of the numerical light-cone we see that if x to +At and Xt -At are 
inside it then 



X ta +At = X a (At) + D_| 



Xto- At = x a (~At) + D- 



with 



D+,D- e 



cAt cAt 



(111.47) 



(111.48) 



The coefficients A and B will then be given by 



1G 



B = -0c 



D A 



D. 



2 (At) 



A = Tc 2 



D A 



D. 



(Aty 



which in turn means 



0- 



!)■ 



i i 

s s 



s (cAt) ' s (cAt) 



(111.49) 



(111.50) 



From this it is easy to see that condition III.43J is indeed satisfied, that is, the three points do form a proper causal 
molecule. Moreover, the absolute value of the acceleration in the new local coordinates will be bounded, and even 
though this doesn't affect the stability of the finite difference scheme, it does improve its accuracy. 

In order to be able to form proper causal molecules everywhere, we must guarantee that two logically distinct 
conditions hold. If we call the central point at to as the "parent" and the points at to ± At the "children" , then every 
parent must have two children and every child must have a parent: 

1. Every parent must have two children. There must always be at least one grid point in the upper and lower time 
levels inside the numerical light-cone of any given point in the middle time level. This can be guaranteed if we 
ask for the distance between grid points to be smaller than the spread of the smallest light-cone, that is, 



Ax < 



2c At 
max (s) 



(111.51) 



which implies 



2 p min ( - | > 1 



(111.52) 



(Without loss of generality, we assume in this section that At and hence p are positive.) 

2. Every child must have a parent. All the grid points in the upper and lower time levels must be inside the 
numerical light-cone of at least one point in the middle time level. This requires that the distance between the 
axes of the numerical light-cones of two consecutive grid points must be smaller than the spread of the minimum 
light-cone. 

Let us therefore consider two consecutive grid points x\ and x 2 — x\ + Ax . The distance between the axis of 
their light-cones at the next time level is given by: 



(x a ) 2 (At) - (x a ) 1 (At) 



(111.53) 



Using the definition of x a we find that: 



d+ = 



Ax - c At 



0(x 2 )-0(x l ) 



\icAt) 2 



T(x 2 ) - T(xi) 



(111.54) 



In the same way we find that the distance between the axis of the light-cones at the previous time level is: 



Ax + cAt 



0(x 2 )-0(x x ) 



\{cAtf 



t( X2 ) - r(xi) 



(111.55) 



The maximum of these two is: 



d = 



Ax + cAt 



0(x 2 )-P(xi) 



l(cAtf 



t( X2 ) - r(xi) 



(111.56) 



Let us assume now that both and T are continuous functions. Then we can expand them in a Taylor series 
around the point: 
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X\ + x 2 



(111.57) 



We then find that, to second order in Ax : 



Ax + cAt 



ox 



\{cAtf 



dT A 

7T Ax 
ox 



(111.58) 



where the derivatives are evaluated at the point x . 

The condition that the maximum value of this distance should be smaller than the spread of the minimum 
light-cone is now: 

2c At 



max (d) < 



max(s) 

Using now the expression for d , we can rewrite this as: 

8(3 



2 p min ( — I > 



I + p Ax max 



dx 



— p (Ax) max 



dr 

dx 



(111.59) 



(111.60) 



This is the "no orphans" condition, that every child point should have a parent. Since we want this to be true 
for all grid points, it must hold for all x . 



We call equations III. 52 and III. 60 the causal reconnection conditions: when they are satisfied, one is guaranteed 



that causal molecules can be formed everywhere. 

It is clear that if the derivatives of (3 and T are too large, it will be impossible to satisfy the second of the causal 
reconnection conditions, Equation III. 60 (no orphans). This can be avoided if we require that (3 and T change very 
little from one grid point to another: 



Ax max 



80 

dx 



< 1 



Ax max 



8T 

dx 



< l 



(111.61) 



As we mentioned before, if this is not the case our finite difference approximation is unlikely to be good anyway, 
and a more refined grid spacing should be used. 

Another interesting feature of condition III. 60 is the fact that, whenever T is not uniform, there will always be 
a value of p large enough for the condition to be violated. These sets an upper bound on the time-step, which can 
be understood if we examine the effects of a non- uniform acceleration on two adjacent numerical light-cones. If the 
acceleration increases with x , the numerical light-cones will eventually converge and pass through each other at a 
large enough time, even if they were diverging initially. Similarly, if the acceleration decreases with x , the numerical 
light-cones will eventually diverge, even if they intersect each other for a while. Clearly these situations do not arise 
in the exact (differential) case because they would break the causal structure of the solutions. We must therefore 
conclude that the numerical light-cones will not approximate the real light-cones properly when we have a time-step 
large enough for these effects to occur. However, if T is such that III. 61 holds, then the upper limit on At will be 
very large indeed, much larger than the Courant limit, and so large that the accuracy of the integration must be 
breaking down anyway. Moreover, for any given time-step condition III. 60 can always be satisfied for a small enough 
grid spacing Ax . 



3. Numerical overheads of causal reconnection. 

Notice that causal reconnection does not change the fundamental structure of the difference equation, since it does 
not affect the relations between points at the final time-step. Therefore, even with causal reconnection, the equation 
will have the form 



Q, 



«+l 



/'( 



(111.62) 



where /' is a different function, which reflects the fact that causal reconnection identifies different points at time-steps 
j and j — 1 to use to generate the points at the final time-step. Therefore, any algorithms that are used without causal 
reconnection for the solution of this tridiagonal system of equations can be used equally well with causal reconnection. 
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There will, of course, be an overhead associated with the search for causally related grid points. In a one-dimensional 
problem with N grid points, this will require only 0(N) operations, since once the causal molecule of one grid point 
has been constructed, the causal molecule of its neighbor will, by continuity, usually differ by at most one spatial shift 
at any time level. In more than one dimension, the search should still scale linearly with the number of grid points, 
since again by continuity the causal molecule of any point can usually be guessed from that of any of its neighbors. 

We have found that for the one-dimensional wave equation, the implementation of causal reconnection given in 
Appendix B can multiply the computation time by something like a factor of two. But for a more realistic problem, 
such as general relativity, where there are many dependent variables per grid point, the overhead of searching for the 
causal structure will be no different than for the simple wave equation, so it will represent a small percentage of the 
overall computing time. 

D. Numerical examples of causal reconnection. 

As an example of the methods that we have developed in the last sections, we will consider a grid that is oscillating 
in the original coordinates (£, t) . The scale and shift functions are given by: 

s(x,t) = 1, 0(x,t) = A cos(u>t) , (111.63) 

from which we deduce 

V = Auj sin(wt) . (111.64) 

This grid turns out to give a very good illustration of all the properties we have mentioned so far. In the calculations 
we have taken c — \ . 

In Figure Owe show two calculations using the finite difference approximation given by Equation III. 15. In these 



examples we have not taken into consideration the causal structure. The figures show the evolution of a Gaussian 
wave packet that was originally at rest at the center of the grid. In both cases we have taken 10 = 6, and we show 
the situation after 35 time-steps (3.5 periods of oscillation of the grid). 

In the first graph A = 1: the maximum shift equals the speed of the waves, but for all the rest of the time the 
shift is less than 1. At the end of the calculation there is no evidence of any instabilities. In fact, we have integrated 
this for a very large number of time-steps with the same result: no instabilities appear. 

For the second graph we have takenA = 1.1: the maximum shift is now slightly larger than the speed of the waves, 
although even here the grid spends most of its time at speeds less than 1. By the end of the calculation an instability 
has started to form. It exhibits the characteristic feature of finite-difference instabilities, that the shortest wavelengths 
are the most unstable. 

In Figure |l2| we compare the direct, non-causal, approach and the causal approach for a larger shift amplitude. We 
use the same initial Gaussian wave packet as before and take A = 1.3, lu — 6. Here the instability appears very fast 
in the direct approach (after only 25 time-steps). With causal reconnection, however, the calculation remains stable. 
We have carried out the same calculation for many more time-steps, and also for larger values of A (up to A — 15 ), 
and the results are the same: no instabilities develop in the causal approach. 

Therefore causal reconnection of the computational molecules seems to cure all the local instabilities on grids that 
move faster than the waves. We will see that in more than one dimension it will also guarantee local stability for 
rapid shifts, but only after we cure a further instability that arises in operator-splitting methods for small velocities. 

IV. THE MULTI-DIMENSIONAL CASE. 
A. How to design an ADI scheme for a hyperbolic system. 

1. Fully implicit scheme. 

We shall begin our discussion of stable integration schemes in more than one dimension by introducing ADI schemes 
in a way that makes our time-symmetric ADI method emerge naturally, and which makes it clear how to generalize 
it to other hyperbolic systems in a straightforward manner. ADI is basically a device for implementing an im- 
plicit integration scheme in many dimensions without the enormous computational overheads that the direct implicit 
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FIG. 11. Oscillating grid, non-causal approach. 
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FIG. 12. Oscillating grid, causal reconnection. 
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scheme would involve. We begin our discussion, therefore, with an examination of the direct implicit scheme and its 
computational demands. We shall concentrate on two dimensions, but the generalization to more is straightforward. 
The general wave equation (Equation [1.3) in two dimensions is 

d 2 (j) r .,., ,~„2l &± 
dy 2 
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(IV.l) 



c dydt dx dy c 2 dt 2 

In the finite difference approximations to this equation, we will use the notation 

<P := <f>l,i y ; (IV-2) 

that is, we will suppress the spatial indices and write them only when the possibility of confusion arises. The finite 



difference approximations to all the differential operators that appear in Equation IV. I have the same form as in the 
one dimensional case, except for a new term that did not exist before: 

8 X 5y(j) 



(d x d y cj)) J 



-t^xy : 



4 (As) 
where the spatial differences are defined in the same way as before, and the truncation error is: 
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(IV.3) 



(IV.4) 



As we learned to do in the one-dimensional case, we will use only explicit approximations for the first spatial 
derivatives. We can then write our second-order implicit finite difference approximation to Equation [V.I in the form: 
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(IV.5) 



In this equation we have assumed for convenience that the spatial increment is the same in both directions, and we 
have defined the Courant parameter in the same way as before. As in the one dimensional case, to arrive at the last 
expression we have multiplied through by (At) . The truncation error is therefore of order: 
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Equation IV.5 is the most direct finite difference approximation to the original differential equation in 2 dimensions. 
We call it the "fully implicit" scheme. As in the one-dimensional case, it takes the form of a matrix equation 



Q 2 ^' +1 = fW, 4S- 1 ). 



(IV.7) 



However, as is well-known, the numerical solution of this equation is considerably more time-consuming than in the 
one dimensional case, and the computational demands increase very rapidly with the number of dimensions. This 
is due to the fact that, if we have N grid points in each of n spatial directions, the matrix Q n will have N n rows 
and columns. Most importantly, this matrix will not be tridiagonal: it may be possible to arrange that the nearest 
neighbors in, say, the x-direction of any point should occupy adjacent columns, but those in other directions will be 
far away in another part of the matrix. The matrix will still be sparse, but the number of operations involved in 
solving it may be very large indeed, in the worst case involving of order N 3n operations at each time-step. Even if a 
well-designed relaxation method is used, the number of operations will in general increase faster than N n . 

ADI schemes offer a systematic way around this problem, usually affording considerable savings in computational 
effort with, as we will show, no sacrifice in accuracy. However, while the fully implicit scheme may be expected to be 
as stable against grid shifts in n dimensions as in one, this is not true of ADI schemes, and we shall have to be careful 
to design a stable one. 
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2. Designing ADI schemes: how to make the operator factorizable. 



Alternating Direction Implicit (ADI) methods W reduce the numerical work involved in an n-dimensional problem 
by modifying the finite-difference scheme in such a way as to replace the original large sparse matrix Q n by one that 
can be factored into a product of tridiagonal matrices related to Q x for each spatial direction. If wc assume that we 
have the same number N of grid points in all directions, we will have to invert a series of N n ~ x tridiagonal matrices 
of size N x N for each spatial dimension. This means that we will need only 0(nN n ) operations to solve the system. 
We see then that the number of operations for the ADI scheme will scale with the number of grid points in the same 
way as it does for an explicit method. 

The reason that one can contemplate replacing the original operator Q n with a different one is that the fully implicit 
finite-difference equation is only an approximation to the differential equation, so if we modify it by adding extra 
high-order terms that are of the same order as those neglected in the original approximation, the accuracy of the 
scheme will not be affected. If we can then choose these extra terms to change the operator acting on the function at 
the last time level <j>> +1 into a factorizable one, we will have speeded up the solution by a huge amount. 



For our two-dimensional wave equation, the operator acting on (j^ +l is (see Equation IV. 5): 



Q 2 := - 1 + | ((3 X S x + (3 y S y ) 



(PT s 2 x + 



n vv 



(pyy 



}■ 



(IV. 



We want to add high-order terms to this expression to transform it into a product of one-dimensional operators of 
the form of the similar term we had in the one-dimensional case (Equation [11.15 with 82 = ) Pc 



So = Q x &/ 



i-2£fc 



'§ 



(try 



2 v 
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(pyy 



(IV.9) 



Let us define S to be the difference between these operators: 



S := Q'n - Q 2 . 



(IV. 10) 



Then we can rearrange the fully implicit equation (Equation IV. 7) to read 



Q / 2 ^' +1 = /(^,^- 1 )+<S^ +1 . 



(IV. 11) 



Now, this is not directly any help, since although we have the factorizable operator Q' 2 on the left-hand-side, we have 
unknown terms in cfP +l on the right. However, let us consider the following related equation: 



Q!^ +l =f{^^- x )+S<tP. 



(IV.12) 



This equation is in a form that can be solved easily, since the unknown (fp + appears only with the operator Q ' . 
Moreover, since in the hmit At — > we have J — * (j>* +1 , in that limit Equation [V.12 approaches Equation tV.ll, 
which is the original fully implicit equation. This fully implicit approximation i s itsel f only a valid approximation 
to the original differential equation in the same limit, so it follows that Equation IV. 12 also approximates the differ- 
ential eq uation in that limit. (It need not be as good an approximation, of cours e: the error terms of the original 
Equation IV. 11 may be smaller than those introduced by the change to Equation [V.12 . We will address this point 
below.) 

There is nothing unique about changing Stfi^ 1 to S<j>* to make the equation factorizable. One could change S<j>* 
to any combination of terms that limits to iS^" 1 " 1 as At — v 0. The different ADI methods make different choices of 



3 This form of the factorized operator is not unique, there are many different operators that one could choose instead of Q x Q v . 
See for example pi and M. 
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these terms. We shall see that the standard choices produce equations that are very unstable when the grid shifts, 
but that by imposing the simple physical requirement of time-reversibility one gets a uniquely defined ADI scheme 
that is stable and just as accurate as the original fully implicit method. 

is: 



It will be helpful to write out explicitly what the operator S defined in Equation IV. 10 



S = - 



p 2 p 3 9 



P x 5 x {[gyy-{pyf] s 2 y } + [f-iFf] 6 2 x (0 y s v ) 



X -W x ) 2 } S 2 X { [g yy - (n 2 ] 5l) 



(IV. 13) 



It is important to note here that the expression for the operator S includes terms in which the difference operators 
in the x direction act on the functions g vv and [3 y , because these functions in general depend on both x and y . 
Had we defined Q' 2 as Q y Q x we would have had a different expression for S , because, whenever the coefficients in 
the finite difference approximation change with position, the operators Q x and Q v do not commute. 



3. ADI schemes old and new. 



We first cast the two original, and still standard, ADI methods for the wave equation, Lees' first and second 
methods, into the notation we have used above. They serve to illustrate how our approach to ADI methods works 
on a familiar method, and we will subsequently analyze the stability of these schemes. Then we will introduce the 
scheme that will turn out to be stable for shifting grids, the time-symmetric scheme. 

a. Lees' first method. The most straightforward approach is that introduced by Lees in 1962 ( M, H) for the case 
of the ordinary wave equation on a fixe d grid. It is convenient to describe it in terms of the extra terms that one adds 
to the left-hand-side of Equation IV. 5 to produce a factorizable equation: 



S 



W+ 1 



y'-^ 



This effectively produces the equation 



Q2 



/j+l 



/(^',^- 1 ) + l s^- 



(IV. 14) 



(IV.15) 



This is a simple change from Equation [V.12 . 

It is clear that as At — > the extra term ( [V.14 ) will vanish, and we will recover the original diffe rential equation. 
However, we will see below that the extra terms do not vanish as fast as the errors in Equation IV. 5, which are given 
in Equation IV. 6. Lees' first method is only first-order accurate on a shifting grid. It is known that Lees' first method 
is absolutely stable on a static grid. However, it will turn out to be subject to weak but significant instabilities when 
used on shifting grids. 

b. Lees' second method. Another way to modify the equation is to add instead a second-time-difference term. This 
is known as Lees' second method: 



S U j+1 - 2cj) j 



tf-v 



(IV.16) 



Here again we recover the original differential equation in the limit At — > 0. The result is in fact a linear combination 
of Equations ( IV.12J ) and (IV.15). We will see below that this method does not sacrifice accuracy: the introduced 
terms are of the same order as the original truncation error on shifting grids. Moreover, it is absolutely stable on 
a static grid. However, our stability analysis will reveal that this method shows strong instabilities when used on 
shifting grids. 

c. The time- symmetric ADI method. A more general approach would be to separate the operator S into different 
pieces and use a first time difference with some of them, and a second time difference with the rest. We will call this 
a mixed ADI method. One can try many different mixed methods, but there is one that is natural from the point of 
view of the original differential equation. This is the one we call the time-symmetric ADI method. 

The original differential equation, Equation IV. 1 , has, in common with all fundamental physics differential equations, 
the property of time-reversal invariance. In this case, the equation is invariant if we make the replacements 



-t 



and 



if 



-{? 



(IV. 17) 
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The fully implicit difference equation, Equation IV. 5, also has this property, since the approximations used for the time- 
derivatives are centered differences: they do not bias the direction of time. In our case, time-reversal is implemented 
by the exchange of the time-step indices j + 1 and j — 1. However, replacing the fully implicit operator Qi with 
Q' 2 breaks this invariance, because this change modifies the way that tp 3+1 enters the equation without automatically 
modifying the cj)^ 1 terms in a symmetrical way. 

If we look at the definition of the operator S in Equation IV. 13, we see that it is not itself invariant: it contains 



terms both linear and quadratic in [3 l . Therefore, since Lees' first and second schemes both add terms in which S 
operates on an expression with a definite time-symmetry, neither scheme is time-reversal invariant. What we need to 
do is to separate S into parts S e and S Q that are even and odd with respect to (3 l and then to allow them to act, 
respectively, on even and odd extra terms. That is, we add to Equation [IV. 5| the term 



S e 



/j+i 



2<^ + ^'- 1 ) + S 



W+ 1 



tf-" 



(IV.18) 



where the even and odd parts of the operator S are defined by 



S, := - 



f3 x 5 x {f3 y 5. 
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(3 X S X 
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g xx -(P x f 6 2 x (JP6 y ) 



(IV.19) 
(IV.20) 



This effectively ensures that we apply the same modification to the ifp -1 terms as to the <fP +l terms in producing 
a factorizable equation that limits to the fully implicit equation as the time-step goes to zero. We will see that, by 
so preserving the time-symmetry of the original equation, we have also produced a method that is just as accurate 
as the fully implicit method and, perhaps more importantly, is unconditionally locally stable on grids shifting at any 
speed up to the wave speed. 



4- Intermediate values and the implementation of AD I schemes. 



Whichever ADI method we choose to use, we will always produce an equation of the form: 



V_ A7 + 1 - 



= At 



Bdp- 1 



(IV.21) 



where A and B are spatial finite difference operators whose specific form will depend on the method chosen. Looking 
at the definition of Q! 2 in Equation IV. 9 , we see that the last equation can be decomposed into a system of two coupled 
equations in the following way: 



pp 



- — Sy-P 2 - 

i PP° x * e 



(a v r 



(p x y 
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fij+1 



■ J+l 



k* 3 + 1 



= At 



BdP-\ 



(IV.22) 
(IV.23) 



i*j+i 



where the first equation defines the so-called intermediate value 

These two equations give the simplest ADI split of the finite difference approximation. In order to solve the system, 
one first solves the second equation for tp* 3 using values of <fi in the previous time levels. This operation involves 
solving a tridiagonal system of equations for each fixed value of the y-'mdex. One then solves for cj>> +1 using the first 
equation, again solving only tridiagonal equations. In the general case of an n dimensional problem, this procedure 
will take us to a system of n equations and n — 1 intermediate values. Each equation employs an operator acting 
only in one of the spatial directions. 



It is important to realize that the splitting of Equation IV.21 given above is by no means unique. One may find 



many different splittings of the same equation, and some may prove to be more computationally efficient than the one 
we have given above. However, the differences will only be in the algebra (and in roundoff errors): different splittings 
are only different ways of writing the same ADI scheme. 
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5. Accuracy of ADI methods. 



The different methods of forming A and B will differ in general in their accuracy and stability. To find the accuracy 
of the different ADI schemes on shiftin g grids, we start by considering Lees' first method. In this case we must add 
to the left hand side of Equation IV. 5 the following term: 

p+i _ 0j-i) = _ L p» Sx [p> s y (&+ 1 - ^'- 1 )] 

-,yy _ (pvy ) s 2 (M+ 1 _ aj~ 1 \ 



S\ 



(g xx - (P x f) si 
P A e 2 



[F Sy (^" +1 - «^- 1 )]} 



(? xx - (n 2 ) s 2 x [(gyy - {n 2 ) % (^' +1 - ^- l ) 



(IV.24) 



The order of this term is found by replacing differences with derivatives: 
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(IV.25) 



From this we can see that the principal part of the truncation error int rodu ced by the new terms is of order (At) , 
which is in fact one order less than the original accuracy of Equation IV. 5 . Using Lees' first ADI decomposition 
reduces the accuracy of the original scheme. This is only true, however, when we consider a moving grid. From the 
last expression it is clear that for a fixed grid (j3 x = (3 y = 0) , the truncation error introduced by this method will 
only be of order (At) , as is well known. 

When we do the same analysis for the case of the ADI scheme based on Lees' second method, we find that the 
principal part of the truncation error introduced by the new terms is of order (At) for a shifting grid. The accuracy 
of the original equation is therefore preserved. In principle, one would therefore prefer Lees' second method. However, 
we will see below that the second method is far more unstable than the first when the grid shifts, so its higher accuracy 
is of limited usefulness. 

For the time- symmetric ADI scheme, the principal part of the truncation error introduced by the extra terms is 
again of order (At) . This method is thus as accurate as the fully implicit one. We will see that it is also stable. 



B. Local stability analysis. 



We turn now to the all-important question of the stability of the ADI schemes that we have described in the last 
section. In the same way as in the one dimensional case, we will start by studying the nature of the solutions of the 
differential equation (Equation IV.l), and we will again consider the solutions in a very small region around the point 
(x l ,t) , assuming that the coefficients remain constant in this region. 

Moreover, for simplicity we will assume that we can take the functions g vv and (3 y out of the difference operators 
in the expression for S (Equation IV.20). Again, if these functions change rapidly from one grid point to the next, 
the accuracy of the finite-difference scheme on this grid is probably poor anyway. 



1. Solutions of the differential equation. 
Following the same procedure as in the one-dimensional case, we will look for solutions of the differential equation 



Equation IV.l that take the form: 



2G 



4> (x,t) = e ia{k)t e lk ' s . 



(IV.26) 



where k is the 2-dimensional wave vector. We denote its components by k l , and define the associated covector 
(one-form) components ki by 



ki — 9ij k 



(IV. 27) 



Substituting into Equation IV. 1 and solving for a, we find the dispersion relation: 

a± = c (kjft) ± c [(kjk j ) + i (k^)] 1/2 . (IV.28) 

This is a simple generalization of Equation III.20J . Again we have the analytic boundedness condition 

\e ra + t e la - t \ 2 = 1. (IV.29) 

Therefore, if one of the solutions is growing, the other is dying at the same rate. 



2. Local stability of the difference equations. 

We will now proceed to the local stability analysis of the different numerical approximations. We look for numerical 
solutions to the finite difference approximations of the following form: 



m i {k^n^+kytiy) Ax 



^ m e 



(IV.30) 



where we have used our simplifying assumption that Ay = Aa:. By substituting this equation into any of the finite 
difference approximations, we shall always obtain a quadratic equation for tp of the form: 



A<f> 2 + Bip + C = 0. 



(IV.31) 



The coefficients in this equation will depend on the particular approximation used. We call the two solutions of this 
equation tp± . 

As in the one-dimensional case, we define the numerical amplification measure: 



M Num := max 
k 



(|V+I 2 ,IV>-I 2 ) , 



and we take our local stability condition to be: 



< m Ae 



(IV.32) 



(IV.33) 



One general cons iderati on applies to all difference schemes. It is not difficult to see that the analytic boundedness 
condition E quatio n IV.29| will also hold in the finite difference case when the following condition on the coefficients 
in Equation [V.31 is satisfied: 



\A\ = \C\. 



(IV.34) 



a. Lees' first scheme. Lees showed that his methods are stable for all time-steps if 9 > 1/2, for a static grid. In 
the shifting case, the coefficients of the quadratic equation for Lees' first method are: 
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{ 1 - ip/3 x sin(k x Ax) - p 2 



(P x 



[cos (k x Ax) — 1] [ 



x [l - ipP v sin(k y Ax) - p 2 9 g vv - {(3 V ) 2 [cos{k y Ax) - 1]} , 

- (V {1-6)1 \g xx - (P x f] [cos (k x Ax) - 1] + \g m - {(3 V ) 2 ] [cos (k y Ax) - 1]} 

- 2p 2 (g xy - /3 x (3 y ) sin (k x Ax) sin (k y Ax) 

- ip 2 Ax [T x sin (k x Ax) + T v sin {kyAx)] + 2) , 



(IV.35) 



(IV.36) 
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c 



(h - ipl3 x sm.(k x Ax) - p 2 9 g xx - {[5 x f [cos{k x Ax) - 1]} 

x |l - ip[3 y sin (k v Ax) - p 2 6 g m - {[3 y ) 2 [cos{k y Ax) - 1]} 
+ 2p 2 9 l\g xx - (/T) 2 ] [cos(k x Ax) - 1} + \g yy - (f3 y ) 2 ] [cos(k y Ax) - 1]} 



(IV.37) 



It is clear that these coefficients do not satisfy the boundedness condition Equation IV. 34 



amplification measure from the roots of Equation IV. 31 



We shall test the stability condition given by Equation IV. 33 on Lees' first method by numerically calculating the 

For simplicity, we consider the case where: 
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(IV.38) 



We do not believe this restricts the generality of our conclusions: from our analysis of the one-dimensional case, the 
restriction on T l should not cause a problem, and the particular values of the metric tensor are unlikely to have a 
determining effect on stability. 

The results of our local stability analysis appear in Figure |l3j, where we show the following region of the shift vector 
space: 

[3 x ,[3 y G (0,1.2). 

Since /3 = 1 corresponds to a grid shifting with a speed c , the region considered in the graphs will include grids 

that shift faster than the waves. We have considered 50 x 50 uniformly spaced values of the shift vector inside this 
region. For a given point in the shift vector space, we find the maximum value of the quantity: 



R := 



AL 



Num 



M An 



using 10 x 10 different values of the wave vector k : 

k x , k y G (0,2tt) , 

and, for each wave vector, 100 different values of the Courant parameter p in the interval (0, 10) . 

Having found i? m ax we plot its value on a logarithmic scale. In the graphs, values of log 10 (-R m ax) smaller than 
or equal to (i? m ax < 1) are represented by clear regions, and values larger than 2 (i? max > 100) by the darkest 
regions. It is not difficult to see that the clear regions will correspond to values of the shift vector for which the finite 
difference scheme is locally stable (at least for p G (0, 10) ), and dark regions to values of the shift that give rise to 
instabilities. The darker the region, the more violent the instabilities. Finally, the solid arc in the graphs marks the 
end of the light-cone. 

It is important to note that the presence of a dark region does not mean that for the given value of lpvec/3 the 
scheme will be unstable for all p G (0, 10) , but only that we must expect instabilities for at least some values of the 
p in that interval. 

In the upper graph in Figure O, we show the case where 9 = 1/4 . The finite difference scheme is unstable for at 
least some value of p at all values of the shift vector. This instability becomes much stronger whenever one of the 
components of the shift vector is greater than the speed of the waves. We recall that even in the one-dimensional 
case, the implicit scheme for 9 = 1/4 is only conditionally stable, so the behavior here is no surprise. This figure 
also shows how the introduction of an operator splitting has broken the rotational symmetry of our problem: it is no 
longer the light-cone which is the important feature, but the rectangular region in which the light-cone is inscribed. 

The lower graph corresponds to the case 9 = 1/2 , which is unconditionally stable in the one-dimensional case. In 
two dimensions, the scheme is still locally stable for values of the shift vector along the direction of the coordinate 
axis, just as the 1-D scheme was. However, instabilities appear for speeds in other directions. These instabilities are 
weak compared to those for speeds faster than the wave speed, but their presence will nevertheless be significant, as 
we will show in the examples of numerical integrations that we give below. 

We have looked at larger values of the parameter 9 , but the situation doesn't improve beyond 9 = 1/2. Lees' first 
scheme is therefore not very useful for grid speeds that are not aligned with the coordinate axis. 
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FIG. 13. Stability for a method of Lees' first type. 
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b. Lees' second scheme. We next turn to Lees' second method, for which the coefficients of the quadratic equation 
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(IV.39) 
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Again, these do not satisfy the boundedness condition Equation [V.34 . In Figure [lj we again portray the cases for 
6 = 1/4 and 9 = 1/2 . The situation is even worse than before: the instabilities in Lees' second scheme grow faster 
than for the first scheme, and even for 9 > 1/2 there is only a very small region of stability just around the origin. 
Clearly, this scheme will not be practical for any moving grid. 

c. The time- symmetric scheme. We have found that both standard A PI met hods become unstable when the 
reference frame is moving. Neither satisfies the symmetry condition Equation IV. 34 . Now w e look in the same way at 
the time-symmetric scheme. The fact that this scheme does indeed satisfy Equation IV. 34 can readily be seen from 
the form that the coefficients of the quadratic equation take in this case: 
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(IV.42) 
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C = h + ipj3 x sin (k x Ax) - p 2 9 g xx - (f3 x ) 2 [cos (k x Ax) - 1]} 

x jl + ipj3 y sm(k y Ax) - p 2 9 g yy - (f3 y ) 2 [cos(k y Ax) - l]\ . 

Figure Q~5^ shows the local stability analysis for this scheme, where again we show what happens for 9 

= 1/2. 



(IV.43) 

(IV.44) 
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FIG. 14. Stability for a method of Lees' second type. 
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ADI: TIME -SYMMETRIC 
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FIG. 15. Stability for the time-symmetric scheme. 
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For the first case, the situation is no better than before: the scheme is unstable for practically every value of the shift 
vector. However, when we set 9 = 1/2, the value that gave absolute stability in the one-dimensional case, the scheme 
becomes locally stable for every value of the shift vector inside the rectangular region that inscribes the light-cone. We 
find that this stability is maintained for larger values of 9 . 

What we see here is effectively a 'light-cone stability condition', except for the fact that instead of a cone we now 
have a rectangle, as a consequence of the fact that the ADI splitting breaks the rotational symmetry of the problem. 
In the general case, this local stability condition can be expressed in the following way: 

g ii _ (^) 2 > o (no sum), (IV.45) 

where i may refer to any spatial direction. 

The time-symmetric scheme has also the important property that in the stable region the numerical solutions will 
always be non-dissipative (at least in the non-accelerating case), that is: 

AfNum 

max ' 



M Ana I V M A 



This can be easily proved from the fact that this scheme satisfies condition [V.34 : If the scheme has a dissipative so- 



lution with AfNum < Af Ana , then condition I V.34 together with the analytic boundedness condition (Equation IV. 2S ) 
implies that it must also have an unstable solution with AfNum > MAna ■ 

With the time-symmetric scheme we have then found what we are looking for: an absolutely locally stable, second- 
order accurate ADI decomposition for the finite difference approximation to the wave equation in a reference frame 
moving at any speed up to the wave speed in any direction. This scheme can be easily generalized to any number of 
spatial dimensions, as can be seen from its definition in the last section. 

The restriction to frames moving slower than the wave speed is expected, of course. To remove it, we now define 
causal reconnection for the 2-dimensional case, using the time-symmetric ADI scheme as our starting point. 

C. Causal reconnection in 2 dimensions. 

In the last section we found that the time-symmetric ADI scheme had stability properties superior to both the 
schemes of Less' first and Lees' second types. However, even for the time-symmetric scheme, instabilities appear as 
soon as one of the components of the grid speed becomes larger than the speed of the waves. In the one dimensional 
case, we saw that this instability could be avoided if we used a computational molecule based on the causal structure 
of the wave equation, and not in the motion of the individual grid points. We now want to generalize this approach 
to the two dimensional case. 

We will again look for a computational molecule that guarantees that the light-cone is properly represented in the 
immediate vicinity of the central point. For the moment we will assume that we have already found the points that 
form such a molecule. We then introduce the local coordinate system {x 1 ,t'} adapted to the causal molecule as a 
direct generalization of the one dimensional case: 

**' := (**-<)-■**(*), 

(IV.46) 

t' := t, 
where {x\ , to} are the coordinates of the central point of the molecule, and where 

pi {t) = A i it -to) + Bi {f _ t(j) + x i Q ] (iy 47) 

A i ._ { x\ 0+ At-teto+x\ -At \ ^ (Iy 4g) 

B i . = f xl 0+ At-^\ -At \ (Iy 49) 

In the new coordinate system, the wave equation has the same form as before, except for the substitutions: 
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A 1 (t - t ) + B l 



V 



A 1 

r*- — 

„2 



(IV. 50) 



Since in the finite difference approximation the coefficients should be evaluated at the center of the molecule, the 
above expressions will reduce to: 



f? 



I? 



B' 
c 



A' 



(IV.51) 



We know that the original finite difference approximation was locally stable as long as Equation |IV.45| was satisfied. 
This implies that the approach based on the reconnected molecule will be stable if 



!J 



0' 



B' 
c 



> 



Vt. 



(IV.52) 



As in the one-dimensional case, we will say that the given three points form a proper causal molecule if the last 
condition is satisfied. 

We will now define a generalization to two dimensions of the concept of effective numerical light-cone. We do this 
by defining first the axis of this numerical light-cone as the line: 



<{t) 



u t 



(Pc (t - t ) + i Pc 2 (t - t ) 2 



(IV.53) 



The numerical light-cone will then be defined by taking at each time the region covered by a rectangle that is centered 
at the axis and has sides: 



2c At (g") 



1/2 



(IV. 54) 



With this definition, the numerical "light-cone" is really a prism and not a cone. It is not difficult to prove now 
that if the points x\ o _ At and x\ 0+At are inside the numerical light-cone of x\ , then we will have: 



/?' + — g -tfT Ag") , (iv.55) 
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(IV.56) 



This means that the three points do form a proper causal molecule, and also that the acceleration in the new local 
coordinates will be bounded. 

As in the one-dimensional case, we can guarantee that proper causal molecules can be formed everywhere if we ask 
for two conditions: 

1. Every 'parent has at least two children. There must always be at least one grid point in the upper and lower 
time levels inside the numerical light-cones of all points in the middle time level. It is easy to see that this will 
require: 



2pmin(<?") 1/2 > 1 



Vi. 



(IV.57) 



2. Every child has a parent. All the grid points in the upper and lower time levels must be inside the numerical 
light-cone of at least one point in the middle time level. We guarantee this by asking that the light-cones of the 
points in the middle time level should cover completely the upper and lower time levels, in other words that the 
union of the intersections of these light-cones with both the upper and lower levels should be the entire grid. 

Let us consider a square of nearest neighbours in the middle time level. We want their numerical light-cones to 
cover the whole quadrilateral area defined by the points where the axis of those light-cones intersect the adjacent 
time levels. A sufficient condition for this to happen is to ask for the sides of this quadrilateral to be smaller 
than the spread of the smallest light-cone divided by y/2 (this factor arises from the fact that the diagonals of 
a square are v2 times larger than its sides). 

Following now the same procedure as before, we can show that this condition takes the form: 
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—= pmin(g") 1/2 > < 1 + 2 max (d n , ^22) 



max(d 2 i,d 2 



1/2 



max(dn,di 2 ) 



Vt. 



(IV.58) 



where the quantities djk are defined as: 



*jk- 



Ax 



dp 



dx k 



(Ax) 



d 2 p 



dx 1 dx 2 



(Ax) 2 (ar* 
2 I da? 



(IV.59) 



and the maximum should be taken over all values of x l . As in the one-dimensional case, the last condition is 
valid only to second order in Ax . 



Conditions IV. 57 and IV.58| are the causal reconnection conditions in the two dimensional case. They will guarantee 
that proper causal molecules can always be formed. 



D. Numerical examples. 

To test the finite difference methods that we have developed, we will consider two different situations: a grid moving 
with a uniform speed, and a grid rotating with constant angular velocity. 

1. Uniformly moving grid. 

We will first study the case of the grid moving with uniform speed, in order to show the advantages of the time- 
symmetric scheme. If the grid is moving with velocity v — (v x , v y ) , it is not difficult to see that: 

(IV.60) 
(IV.61) 
(IV.62) 

Using these values for the coefficients, we have studied the numerical solution to the wave equation for a number of 
examples, comparing the three different ADI methods developed earlier. The first set of graphs (Figure Iq) show the 
result of one such calculation for a scheme of Lees' first type. In the graphs we show the grid region [(0, 10) x (0, 10)] , 
and we calculate the evolution of a Gaussian wave packet originally at rest at the point (7, 7) . For simplicity, we have 
imposed reflecting boundaries. We have taken a time-step such that p = 1, which means that we are well beyond the 
Courant limit {j The evolution is followed using a grid with a speed given by: 



9ij — °ij 




F = v'/c, 


and 


pa; _ py _ 


0. 



V = ( 2'2 } ' 



\v\ = 0.707 < 1, 



where we have taken c = 1 . 

We see how an instability is beginning to grow even though the grid is moving slower than the wave speed. This 
is precisely in accordance with the results of our local stability analysis. This instability grows slowly, as expected. 
Nevertheless, it is clear that its presence is unacceptable in a calculation of any duration. If we use a method of 
Lee's second type, the instability takes longer to develop, but once it appears it grows very fast, much faster that 
with Lees' first method. The fact that the instabilities in general take longer to appear with Lees' second method 
can be traced to the particular wave modes that are involved. As we can see in the graphs, the instabilities in Lees' 
first method are associated with relatively long wavelengths (several grid points), and since these modes are already 
present in the initial data, they start growing right away. In Lees' second method, however, the instabilities turn out 



In a n dimensional problem, the Courant limit for the stability of an explicit scheme is p = 1/yfn. 
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FIG. 16. Uniform shift vector: Lees' first scheme. 
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to be associated with very short wavelengths (one or two grid points), which do not contribute significantly to the 
initial data. These means that, even though the instabilities are more violent with this method, it will take a long 
time for the unstable modes to grow to the scale of the real solution. 

In the next set of graphs (Figure IVm we have applied the time-symmetric scheme to the same problem. The 
instability has completely disappeared. This is again in agreement with our previous conclusions, and shows the 
superiority of this method. 

We have performed similar calculations for many different values of the grid speed and we have found the essentially 
similar results. The time-symmetric scheme remains stable as long as the grid moves slower than the waves, while 
the other schemes present instabilities for quite small grid velocities. 

2. Rotating grid. 

In order to show the advantage of causal reconnection of the computational molecules when the grid moves very 
fast, we will consider now an example with a grid rotating with the constant angular velocity w . It is not difficult 
to show that: 

9ij = 5 tj (IV.63) 

[3* = --y, P> = -x, (IV.64) 

c c 

v x = -^ x , ry = -^ y . (iv.65) 

c c 

To test for local stability when using causal reconnection, we take 

u> = 0.25 

in units in which c = 1 and the grid extends over the range [(—5, 5) x (—5, 5)] . This means that the centers of the 
edges of the grid will be moving faster than the wave speed, with a linear velocity of 1.25, while the corners will be 
going even faster. 

The next graphs show the results of a time-symmetric calculation, first using a "direct" calculation (fixed compu- 
tational molecule) and second using causal reconnection. Again we use a time-step such that p = 1. In Figure ttq, we 
show the evolution of a Gaussian wave packet originally at rest at the center of the grid, using the direct approach. 
We see how after 32 time-steps an instability has appeared close to the boundaries. Only five time-steps later, this 
instability has grown so large that the original wave is no longer visible (the scale is automatically adjusted to display 
the largest value of the function). 

Figure H9 shows the same calculation using causal reconnection. The instability is not present. In fact, we have 
done the same calculation with much larger values of the angular velocity (up to w = 3.0 , where the edge is travelling 
at 15 times the wave speed), and the scheme remains locally stable. 

These examples demonstrate dramatically that time-symmetric ADI can be married with causal reconnection, and 
that together the two techniques provide a robust difference approximation to the wave equation on a moving grid. 
These methods are stable, offer all the computational advantages of ADI schemes, and remain second order accurate 
in Aa; and At . 

A comment on how to enforce causal reconnection at the boundaries seems in order here. In all our examples we 
have taken the practical approach of setting the value of the wave function to zero whenever a complete causal molecule 
cannot be formed. This can happen not only at the boundaries, but also at inner points close to the boundaries for 
large enough grid speeds. The philosophy behind this approach is simple: If the causal molecule is incomplete, then 
we would need information from outside the grid to evolve the wave function at that place. If we impose the condition 
that no information can come from the outside, then we must take the value of the wave function as zero at that 
point. This requirement can be relaxed somewhat be using an outgoing wave boundary condition whenever we can 
still find a causally related point in the previous time level, but not before that. At places where one can't even find 
a causally related point in the previous time level, the only legitimate thing one can do is to set the value of the wave 
function to zero. 
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FIG. 17. Uniform shift vector: time-symmetric scheme. 
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FIG. 18. Rotating grid: non-causal approach. 
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FIG. 19. Rotating grid: causal reconnection. 
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V. CONCLUSIONS. 

The wave equation we have studied here is a prototype for more complex equations of mathematical physics, such as 
the Einstein field equations. In fact, many hyperbolic systems in mathematical physics can be formulated in terms of 
the wave operator. One would expect the instabilities we have found here to be generic: any numerical approximation 
to a hyperbolic system on a shifting grid should exhibit them. 

Only experience will show us just how well our cures for these generic instabilities transfer to more interesting 
equations. However, the instabilities we have described here are cured by the application of two clear physical 
principles, causality and time-reflection invariance. It seems clear that it would be asking for trouble not to incorporate 
these principles into the design of algorithms for the numerical integration of any fundamental physical equations. 

We have, of course, studied in detail only one second-order differential equation in one and two dimensions. The 
restriction to two dimensions is not important. The physical principles involved do not depend on the number of 
dimensions, and the savings obtained by using an ADI scheme instead of a fully-implicit formulation increase rapidly 
with the number of dimensions. In many physical systems, it is advantageous to formulate the equations of the theory 
as first-order differential equations. This is true in hydrodynamics and in many studies of general relativity. The 
general principles of causality and time-reflection invariance extend in a simple way to such systems. Time-symmetric 
ADI should prove relatively straightforward to apply to more complicated systems of equations, provided the original 
differential equations embody time symmetry. 

Causality may be less straightforward in nonlinear equations, where the structure of the characteristic cone will 
depend on the solution, and so the exact causal relationships between time-levels cannot be decided independently of 
solving the equations. However, causal reconnection is implemented via an inequality: one requires that grid points 
should be within the characteristic cones of their relatives at the previous time-step. In most cases, one would hope 
that the inequality can be assured simply by extrapolation from the behavior of the characteristic cones on the known 
time-steps. 

An important area for the application of the techniques we have developed here would be numerical fluid dynamics, 
where the study of wave phenomena in supersonic flows is a natural place to expect causality problems, and where 
the interest in three-dimensional problems makes ADI essential in many cases. 

In some restricted situations, it may be straightforward to apply these techniques; for example, a neutron star moving 
supersonically through a grid in general relativity will, if treated in the standard way, use acausal computational 
molecules. Using causal reconnection, adapted to the characteristics of the fluid problem, should prevent instabilities 
of the type we have found here. 

But the application of our techniques to more general problems in fluid dynamics will not be automatic. Causal 
reconnection will have to be generalized to deal with hydrodynamic shocks. At a shock, the regular causal structure of 
the fluid breaks down. This does not mean that causal reconnection cannot be implemented there. On the contrary, 
the fact that a causal algorithm is constantly mapping the structure of the characteristics means that it can be 
programmed automatically to locate and to identify shocks. 

The idea of correctly representing the causal structure of the original differential equation is not new, existing 
methods for handling shocks and related transport problems, such as upwind differencing ||] and Godunov methods 
JlO| |, are already based on the local structure of the characteristics of the fluid. These ideas have also been introduced 
in the numerical study of steady supersonic flows, where the direction of flow behaves like a time coordinate and the 
equations become hyperbolic. Integration methods have been developed that use retarded differences in the upstream 
direction to maintain stability [|ll| , pg] . All these methods differ from causal reconnection in the fact that they keep 
using only the nearest neighbours to build the computational molecules. We have recently become aware, however, 
of a paper by E. Seidel and Wai-Mo Suen that introduces an idea they call "causal differencing" , that is very similar 
to our causal reconnection p3[ . 

Fluid dynamics also presents special challenges to time-symmetric ADI. The usual equations of inviscid gas dynamics 
are time-symmetric, but the presence of viscosity or shocks introduces a fundamental irreversibility into the problem. 
We hope to treat the fluid dynamic problem in a future paper. 

We are confident, however, that the present techniques will generalize easily to problems in numerical general 
relativity, such as that of the motion of black holes through fixed grids. Causal reconnection should allow the equations 
to remain stable and causal. Moreover, the computational advantages offered by ADI schemes, of permitting stable 
large time-steps (provided the physical situation allows such steps to remain accurate) while avoiding time-consuming 
sparse-matrix solutions, can be obtained without sacrificing accuracy or stability. It is hard now to imagine any 
situation in numerical integrations of the vacuum field equations of general relativity where one would use implicit 
methods without employing time-symmetric ADI. 



41 



APPENDIX A. DERIVATION OF WAVE EQUATION ON A SHIFTING GRID. 



In this appendix we will sketch the derivation of Equation [1.3 by making use of elegant tensorial techniques. There 
are many alternative approaches, of course, and a reader unfamiliar with tensors can obtain the same result in a 
straightforward, but rather long, way simply by making the following general change of variables in the original wave 
equation from the physical (inertial) coordinates {£ M } to the computational coordinates {x a }: 



x % (t 



V«M 



€° : 



(A.l) 



with the associated change of derivatives 
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dx k d 
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dx k d d 
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d£,° dx k ' dx° 



(A.2) 



The functions that we have identified as the shift vector j3 l , the spatial metric gij and the T l coefficients in 
Equations ([1.4), (II. 5) and ([1.8) come out as part of the algebra. A reader who wants an introduction to the use of 
tensors in mathematical physics is invited to consult reference jl4|. 

We will start from the expression of the wave equation in a general coordinate system: 



□ 2 = (-r t^u = o, 



(A.3) 



where the semicolon stands for covariant derivative. Using the explicit expression for the covariant derivatives, the 
last equation takes the form: 



•jW . 
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dx^ dx u 



r A — — = o 

dx* ' 



where the coefficients T x are defined in terms of the Christoffel symbols as: 

pA uv pA 

Our first task is, then, to find the inverse metric 7 Miy . The metric in the new coordinates is given by 



1^ 



dx^ dx v 



Va/3, 



with rj a p the Minkowski metric tensor: 



(A.4) 
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We now note that for a line of constant {£'} we have 
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Using now the definition of the shift vector (Equation [1.4 ) and writing x° = ct we find 



dx° dxJ 



This is an important relation, and we will use it to rewrite the metric coefficients given by Equation A.6 
For the mixed components in space and time of 7 Mt/ we find: 



(A.7) 
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(A.9) 



42 



^ d? d? 

70* - "^^0- l^l^WrO 



1=1 



= t?Y. 



I p,rf 



1=1 



dx % dxi 



and finally: 



7oi = Ho = Hjl3 J = gi]l3 J 



(A.10) 



In a similar way we can find the coefficient 700 : 
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and from this we obtain: 
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(A.11) 



We will adopt the convention that the indices of the shift vector can be raised and lowered by using only the spatial 
metric: 
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(A.12) 



where g tJ are the coefficients of the inverse of the spatial metric matrix gi j 
The coefficients of 7 M „ can now be written as: 
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Using the last expression it is not difficult to see that the coefficients of the inverse metric 7^" will be given by: 

-1 k 
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(A.14) 



J (g Jk -0 j k ) 



Having found 7^", we will now look for an expression for the coefficients T l . Since the original coordinates {£"} 
define an incrtial reference frame, the Christoffcl symbols can be expressed in terms of their transformation to the 
general coordinates: 



A _ dx x d 2 i a 
» v ~ d^ dxt* dx" 



From the last expression it is easy to see that: 



which in turn means: 
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On the other hand, from the general expression for the Christoffcl symbols: 
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it is not difficult to show that: 



■^{a^ + BT^-'W] 



Using the previous results, we can finally rewrite Equation A. 4 in the following way: 

d 2 (J3 2/3* d 2 (j) 



(g ik ~f3 l (3 k ) 



dx l dx k 



c dx l dt 



T id±_} L cP± = Q 



dx i 



dt 2 



(A.19) 



(A.20) 



This is the final form of the wave equation in the coordinate system adapted to the motion of the grid. 



APPENDIX B. IMPLEMENTATION OF CAUSAL RECONNECTION. 



In this appendix, we discuss one algorithm that determines the positions of the points that form the causal compu- 
tational molecules. We will consider the case of an arbitrary number of spatial dimensions n . The particular cases 
of one and two dimensions can then be found in a straightforward way. 

There are many different ways of finding the closest causally connected grid points. For example, if it is possible 
to find the transformation of coordinates that takes us back to the original inertial reference frame 



e = e{x\t) , 



(B.l) 



then we could use the fact that in that reference frame the causal structure is particularly simple: we would simply 
select those points in the different time levels that have the closest values of {£'} . This method, however, will only 
be useful in a few special cases. Indeed, in the general case it may prove almost impossible to find the functional 
relation (B.l). 

With this in mind, we have developed a method that can be applied in the general case. Let us then assume that 
we are given a point in the last time level {to + At} with position x\ , At . Our aim is to find points x\ and x\ _ At 
in the two previous time levels in such a way as to guarantee that a proper causal molecule will be formed. We have 
already seen that this will happen if both x\ _ At and x\ +At are inside the numerical light-cone of x\ . 

Our algorithm to find the proper causal molecules linking the grids at time-steps to + At, to, and t — At assumes 
that the causal reconnection condition holds. (If it doesn't, then remedial action, changing At or Ax 1 , is required.) 
Our procedure is the following. 

1. Choose some point x\ +At . The center of its causal molecule will be that grid point y\ Q for which the following 
function reaches a minimum: 
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(B.2) 



This minimum can easily be found by standard multi-dimensional search techniques. Once the minimum is found, 
we have our best approximation to the exact inverse coordinate transformation: the point x % is approximately 
at the same spatial location as y l in the original inertial frame. 

2. Once we have found the appropriate y\ , we look in the third time level {to — At} for the completion of the 



causal molecule, the point z\ _ At that minimizes the function: 



N ( r i ii 

h (**) := E { [yl + p (vL . *) At + 2 r (H ' *) ( At f\ - zl } 



(B.3) 



This is much easier than the previous step because we already have the point yl , so we don't have to calculate 
again the functions (3 and T . In fact, minimizing ji is equivalent to finding the point z l in the third time 
level that is closest to the axis of the light cone of y\ Q in the original inertial frame. We do not have to go 
through the grid again to find this point. 
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3. So far we have constructed only one molecule. One needs to repeat the above steps for all points at time-step 
to + At, but of course the best guess for a causal molecule for any grid point is simply to translate the molecule 
found for its neighbor. This will occasionally fail, but only by one grid point. So the minimization steps will 
require a computing effort that is only proportional to the number of grid points. 

It is not difficult to prove that, whenever the causal reconnection conditions hold, this algorithm does indeed produce 
a proper causal molecule. For reasons indicated above, the computational effort is proportional to the number of grid 
points. In complex problems, such as general relativity, this is likely to be a very small overhead. 
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